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ff^ I Abstract. This paper compares the results of two SBBN codes developped 

independently by different teams of physicists. These two codes have significant 

^ . differences that lead to a discrepency between their final mass fractions of ^He 

I of 0.003. This paper shows that the mass fractions of each code had different 

' orders of convergence, and how the number of timesteps affects the accuracy of 

\^ ■ the mass fractions. At the end, the paper shows how to modify both codes so 

', that their ^He mass fractions agree to around 0.0001. 

OO . 

a\ ■ 

^ ' 1. Introduction 

Oh' 

^ , In 1993 I started working on a Big Bang Nucleosynthesis code once used by UT- Austin 
^ I Professors Tony Rothman and Richard Matzner ( hereafter known as the Texas code ) 
^ ■ for some papers between 1982 and 1984.0 Two years later Mr. David Thomas e-mailed 
^ ■ me a different BBN code ( hereafter known as the Thomas et al code ) which he used 
• 1— I . with Profs. David Schramm, Keith Olive, and Mr. Brian Fields for a series of papers 
^ i starting in 1993.0 But these two different codes calculated different final mass fractions 
I of their isotopes, as shown in Table (1). In particular, the values of Xi^e disagreed by 
about 0.003. So I went about determining the differences between the codes that could 
explain this discrepency. 

A Big Bang Nucleosythesis code starts with free neutrons and protons plus photons, 
electrons and the three neutrino species, all at a high value of the electron/photon 
temperature T. Then in timesteps of value At the code calculates the current abundances 
Y{i) of each isotope as determined by the reactions between the isotopes and other 
particles, as well as the temperatures and energy densities of the non-baryonic particles, 
which affect the reaction rates. The Thomas et al code featured some details like the 
Coulomb factor correction to the neutron-proton conversion rates that I had to add to 
the Texas code. Also the codes had several differences in calculating various thermal 
quantities.0 The Thomas et al code would for instance calculate T via the Runge- 

^The papers dealt with anisotropic models, but the code used for them had the option to turn the 
anisotropies off. Also I expanded this code's original number of 9 isotopes to the other code's 68. 

^These papers lead up to a 1994 paper on an inhomogenous model, which could be compared to the 
standard model of the Thomas et al code. 

■^Kawano ( 1992 ) describes the details of a code similar to the Thomas et al code. 
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Table (1): Mass Fractions X from each code. The results are for 7710 = 3.0, but the codes 
have the same discrepency between the values of Xi^g for all interesting values of r/io- 



Kutta method and then plug in T to equations for energy densities and pressures. The 
Texas code would use R-K to calculate density pe+7 of the electrons and photons and 
then determine T from there. I wondered which differences caused the disagreements 
between the codes. This paper details two differences in the way each code goes about 
its evolution that made the accuracy of each code's results questionable. 

2. Convergence of the codes 

At each timestep n the BBN codes use the Second-Order Runge-Kutta method to cal- 
culate the isotope abundances 



Wagoner ( 1969, p. 253 ) lists the complicated equation for Y{i), which depends on the 
other abundances and on the reaction rates that destroy ( [ij] ) and create ( [kl] ) isotope 
i. But at the highest temperatures the total destruction rate of each isotope is nearly 
equal to its creation rate. So the codes calculate Y{i) using an equation linearized in 
terms of y^+i(i) 



dt 



dt 



At 



(1) 



-^it,Yn) = A,(i;,[zj](r„),[M](r„))xi;+i(j) 

where the matrix Aij is written out in Wagoner (1969, p. 294 ). The codes pair off the 
right hand side of this equation with an implicit expression of Y{i) leading to a matrix 
equation for 



dY(i) 

Yn{i) = r„+iW-^^[t + At„,F„+i(i)]At„ + --- 

= Qi,{Yn,mT^),[kl]{Tn)) xYn+,{j) (3) 

( Qij = Ijj ~ A'ijAtn ) The codes solve for l^n+i(0 and plug it back in to get Y{i). 

This R-K method should produce mass fractions that have second-order convergenge 
as the timestep values get smaller. I tested the convergences of mass fraction XiHe the 
Thomas et al code and the Texas code by starting with an array of timestep values At„ 
that naturally arise from a run, and then dividing those values by two, four, and so on. 
The Thomas el al code turned out to have the expected convergence: 



XiHe = 0.2392602277 - (0.2474541681 x 10-^)^^ 

where d equals one over the number of divisions of each timestep. But the Texas code 
had only first order convergence, and so needed to be fixed. 



X^He = 0.2394509836 + (0.3832102957 X 10"^)d 

I traced each code to find out exactly how in step # n each went about calculating 
Ynj^i{i). In the first Runge-Kutta step the Thomas et al code used the previous F„_i(i) 
on the left hand side of matrix equation (3), and Yn{i) in matrix Qij. 



Yn-i{i) = Q^J{YnM{Tn)Akl]{Tn))^Yna{3) 

'^^%,Yr,{i)) = Yna{i)-Yn-,{i) 



dt ' ' At 



n-l 



I'm calling the equation solution Yna{i) for the first step ( and then F„^(i) for the 
second step ). Then this code used an expression for Y{i) that clearly matched up 
with Equation (2) for Y{i)[t,Yn{i)], with Yn-i{i) in the numerator and Atn-i in the 
denominator. So Yna{i) appears to be an approximation of Y^. The Thomas et al code 
then uses Y{i)[t, Yn{i)] and a new timestep Atn to get the value Yn{i) 

Ui) = y„(i) + ^^^^^^]rf^^At„ 

an approximation of Yn+i{i) to be used in the matrix of the second Runge-Kutta step. 



y„(!) = 0., >;..|y](r„),[H](T„) xr„^(j) 



dt ' ' ~»-"^^v"// At, 

I ^Y^t) - F„_i(z) ^ M^)-y„(^) j 



y„+i(i) = Ynii) + ^ 
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Atn-i At, 



And this codes used in its matrix equation. So one could say that Y^i^ was an 

approximation of and this second time derivative corresponded to Y{i)[t + 

Atn,Yn+i{i)] according to Equation (2). Plugging them into Equation (1) we get Equa- 
tion (4) for our actual value of Yn+i{i). Note the Atn as a factor and the Atn-i in one 
of the denominators. 

The first R-K step of the Texas code, though, has Yn{i) on the left-hand side, and 
Atn in the denominator for Y{i).: 

Yn{i} = Qij{yn,miTn),m{T„))xYn^{j) 

dt ^' " Atn 

That would imply that this code's Ynaii) is an approximation of Yn+i{i), and that 
Y{i){t, Yn) has been determined explicitly instead of implicitly. The second R-K step, 

Yn(i) = Q«(V;,lijl(T„+i),[H](T„+i)) xy„^(j) 

dt ^ ' ' Atn 

implies that Ynfs is also an approximation of Yn+i, and when the Texas code plugs in our 
expressions for the F(i)'s.: 

YnMi) = yn{i) + l 

= l{Ynaii)+Ynp{i)) (5) 

the Atn factor cancels out the At„'s in the denominators. So the final Yn+i{i) is the 
average of the two first-order Euler approximations Yna and Yn/3, calculated from solving 
the matrix equations alone. 

So I modified the Texas code to follow the Thomas et al code's more consistent scheme 
of calculating The Texas code's Xah^ then acquired the second-order convergence 

it should've had.: 

X^He = 0.2394508633 - (0.1509902093 x 10-^)d^ 

For d = 1, X4He fell from 0.23970 to 0.23929. Still far from the other code's 0.23690. 
But X4He in the Texas code converged to the same value as it did with first-order. 

3. The value of Ymin 

Y{i) of the Texas code could go all the way down to 0.0. But the Thomas et al code 
put a lower limit Ymin — 10~^^. That code used Y^in in its equation for Atn 
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Atr 



+ 



Atr. 



Atr. 



Tin GK 

Figure la; Timesteps At^ of the Thomas et al code for 
various values of minimum abundance Y^|„. The number of 
timesteps greatly increases with decreating Y^^,^. The linear 
portions of the graphs correspond to when the code uses 
temperature T to determine At„. 
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Figure lb: Mass fraction X.jjg of the Thomas et al code 
as a function of Y^^. The total number of timesteps varies 
from 323 for Y^|„ = lO-^^ (.q 1020 for Y„,„=10-85. 
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The Thomas et al code compared the ratio of T to T to the Y{i)/Y{i) ratios to pick the 
smallest ratio times dk = 0.1 as the value of Atn, so that the code could evolve stably. 
But to the Y{i)/Y{i) ratio the code also put a factor dkfac, designed to prevent the code 
from picking a Y{i) whose value was close to Ymin- But this factor lead to larger values 
of Atn than in the Texas code, and hence much fewer timesteps. At 7710 = 3.0, the 
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Figure 3a: The Thomas et al code it„ now with the 
temperature ratio T/ T lowered by a factor of 0.1. This 
ratio now determines the value of it^ for nearly the 
entire run, regardless of Y^j„ 
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Figure 2b: With this new 0.1 factor, the dependence of 
X.jjg on Yjjjj^ is greatly reduced in the Thomas et al code. 
Graphs very similar to Figures (l — 2 ) can be made from 
the Texas code. 



Thomas et al code would run in only 323 timesteps, whereas the Texas code needed 
2772. Figure (la) graphs At„ for various values of K„jj„. Varying from 10^^^ to 
10^^^ would vary the number of timesteps from around 300 to 1000. And Xa^^ would 
vary from 0.257 to 0.259, as shown in Figure (lb). 

The hnear sections in Figure (la) corresponded to when the code chose T/T to 
determine At„. So I got the idea of putting on T/T a second factor 0.1. Figure (2a) 
shows that this 0.1 factor lowered the T/T ratio to the point that the code would choose 
that ratio nearly all of the time. So in Figure (2a) the number of timesteps varied only 
from 1000 to 1300. And Figure (2b) shows that the 0.1 factor eliminated the influence 
of Y^in on X4He. So the number of timesteps instead of l^rnm itself really determined the 
accuracy of ^4^5. 

The Texas code uses this equation to determine the timestep. 
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No dkfac factor. So the code could run with Ymin set at zero. I modified the Texas code 
to have a dk^ac for a non-zero Ymin- I also wanted to put in a 0.1 factor, but instead of 
T the Texas code used which depended on T^ for most of a run. So I put in a 0.4 
factor instead. And indeed the Texas code exhibited the same behavior as the Thomas 
et al code with these new conditions, though Xahc tended to vary more here than in the 
other code. 

I'd still prefer to not have a non-zero Ymin at all. But the Texas code's old way of 
calculating Ai„ would have that large number of steps in that case. So for Ymin — 0.0 I 
put in both codes the following dtfac factor.: 



f QiC J J, 

dk 

This factor set the number of timesteps at 1320 for both codes. A number on the 
order of 1000 seemed sufficient to get a value of Xi^e that didn't depend very much on 
the timestep number. 

4. Final results 

Table (2) shows the final mass fractions for each code with the modifications put in, for 
the case of 7710 = 3.0. The codes now have close agreement with each other, especially 
for ^He where the mass fractions are within 0.0001 of each other. These results are for 
Ymin = 0.0, but I got similar results for Ymin = 10~^^ as well. 

As for the thermal quantities, those differences seemed very confusing. But I checked 
these thermal quantities and determined them to be nearly equal between the codes after 
I put the changes in. So the number of timesteps and the convergences were the most 
signficant reasons for disagreement between the codes. 
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Table (2): Mass Fractions for Modified Codes Here, Ymin = 0.0 and 7710 = 3.0. But the 
codes still agree for the range of 1 to 10 for 7710 and Ymin up to 10~^^. 



And finally, I took convergence tests of the codes again after all the changes had 
been added. The Thomas et al code obeyed: 

XAHe = 0.2392673182 - (0.3253724834 x 10~^)d^ 
while the Texas code obeyed. 

X^He = 0.2394508628 - (0.1203603394 x 10"^)^^ 

The codes still converged to the same values that they've converged before. So the 
remaining differences between the codes result in these final values disagreeing by 0.0002. 
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